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Abstract 

We present here a semi- analytical solution of the problem of particle acceler- 
ation at non-linear shock waves with a free escape boundary at some location 
upstream. This solution, besides allowing us to determine the spectrum of par- 
ticles accelerated at the shock front, including the shape of the cutoff at some 
maximum momentum, also allows us to determine the spectrum of particles 
escaping the system from upstream. This latter aspect of the problem is crucial 
for establishing a connection between the accelerated particles in astrophysi- 
cal sources, such as supernova remnants, and the cosmic rays observed at the 
Earth. An excellent approximate solution, which leads to a computationally 
fast calculation of the structure of shocks with an arbitrary level of cosmic ray 
modification, is also obtained. 
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1. Introduction 

The importance of the process of particle acceleration in astrophysical shock 
waves for the origin of cosmic rays (CRs) is now generally acknowledged but 
several weak points remain in the theory when one tries to establish a con- 
nection between accelerated particles and cosmic rays observed at Earth. The 
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main problem is related to the difficulties at assessing the role of escaping par- 
ticles: while supernova remnants (SNRs) are often invoked as the main sources 
of Galactic CRs, at least up to the knee, their ability to generate CRs with 
the spectrum observed at Earth is all but proven. If particles are trapped 
in the expanding shell during the Sedov- Taylor phase, adiabatic energy losses 
prevent the release in the interstellar medium of particles with energies in the 
knee region. If SNRs are the main sources of CRs up to the knee, ongoing 
escape of particles from the upstream region is required during the Sedov- 
Taylor phase. The spectrum of these escaping particles is expected to be peaked 
around the maximum momentum that can be reached in the accelerator at a 
given time. The actual spectrum of CRs from an individual SNR is therefore 
the convolution over time of these peaked spectra. Despite the crucial impor- 
tance of escaping particles, their role in cosmic ray modified shocks has re- 
ceived scarce attention so far, with some noticeable ex ce ptions (see for instance 
the w o rk by Bcrczhko, Yclshin fc Ksenofontov ( 1994) 



Ptuskin fc Zirakashvili 



(2005); 



Lee. Kamae fc Ellison 



2008); 



Reville et al 



20081) : 



Caprioli. Blasi. Amato 



20091 )). One of the difficulties from the technical point of view is that it is not 



clear which particles do actually escape the system. While from the mathemat- 
ical point of view, escape can be modeled by requiring the existence of a free 
escape boundary upstream, from the physical point of view the issue remains 
that the position of this boundary is related to poorly understood details of the 
problem, especially the ability of particles to self-generate their own scattering 
centers. The position of the free escape boundary could coincide with a location 
upstream of the shock where particles are no longer able to scatter effectively 
and return to the shock. This however would lead to an anisotropic distribu- 
tion function of the accelerated particles, that can no longer be described by the 
standard diffusio n-convec t ion eq uation. Moreover, wh ile waves ca n be generated 



both resonantly (|Skilline 



19751 ) and non- resonantly (Bglj , 120041) , particles can 



scatter effectively only with resonant waves. This adds to the complexity of the 
problem, in that one might have amplified magnetic fields of large strength but 
on scales which do not imply effective scattering of the highest energy particles. 
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In the absence of a better description of this phenomenon, so far the best 
way to mimic the escape is to impose a reasonable location for a free escape 
boundary and calculate the escape flux as derived from the diffusion approxi- 
mation. Here we present an exact semi-analytical solution of this problem for 
shocks with arbitrary cosmic ray induced modification. We also propose a sim- 
ple approximate expression which turns out to be an excellent approximation to 
the exact solution. The approach presented here allows us to calculate the spec- 
trum of accelerated particles at any location upstream and downstream of the 
shock and the spectrum of escaping particles, in the assumption that a quasi- 
stationary situation is reached at any given time. Clearly within this approach 
the maximum momentum achieved by the particles is not imposed by hands but 
rather obtained self-consistently from the condition of free escape at xq. 

This is not the first attempt in the literature at investigating the problem 
of free escape from a sh ock r egion: the problem was re c ently faced numerically 
by iReville et al.l ([20081 ) and IZirakashvili &: Aharonianl (12010 ), who spec i alized 



their calculation to the case of SNR RX J1713. 7-3946. 



Kang fc Jones! ( 1995 



20061 ) investigated the problem of particle acceleration at a modified shock with 
free escape through a time-dependent finite differences scheme with adaptive 
mesh refinemen t of th e grid. A Monte Carlo technique w as adopted by e.g. 



Jones fc Ellison 



(|199ll ); IVladimirov. Ellison fc Bvkovl (|2006l ) to have a handle 
on the escape flux of particles from the shock. It is worth stressing however that 
these numerical methods require computation times for a given set of parameters 
which range between several hours, for the Monte Carlo technique, and several 



days for the time dependent calculations of iKang fc Jonesl (|1995l 120061 ). These 
times should be compared with a typical computation time of 1-2 minutes (on 
a laptop) required for the semi-analytical method discussed here or pr evious 
versions of it, in which a boundary con dition in momentum was adopted (Blasi, 



2002 


2004: 


Amato & Blasi 


2005 


2006) 



2006). The issue of computation time becomes 



crucial when these calculations are embedded in hydrodynamical codes for the 
evolution of SNRs. 

The paper is organized as follows: in Sj2]we obtain the implicit exact solution 



3 



of the problem with a given free escape boundary condition; the approximate 
solution, along with its comparison with the exact one, is presented in [3] We 
conclude in fQ] 



2. Exact solution 



We start from the stationary, non-relativistic, one dimensional diffusion- 
convection equation for the isotr opic part of th e distribution function of accel- 
erated particles f(x,p) (see e.g. 

df(x,p) d 



Skilline 



1975): 



u(x)- 



dx 



dx 



D(x,p) 



df(x,p) 
dx 



p du(x) df(x,p) 



3 dx 



dp 



+ Q(x,p), (l) 



where D(x,p) is the diffusion coefficient, with arbitrary dependence on both 
position and momentum, Q(x,p) is the injection rate and u(x) is the fluid ve- 
locity in the shock frame. Here, for the sake of clarity, we neglect the velocity of 
the scattering centres, which is typically of order of the Alfven velocity va, with 
respect to the fluid. The generalization to the case of small Alfve nic Mach num 



ber m ay be easily obtained following the procedure discussed in (|Caprioli et al 



20091 sec. 3). We solve this equation with the upstream boundary condition 
f(xo,p) = 0, which mimics the presence of a free-escape boundary placed at 
a distance xo upstream of the shock (placed at x = 0). The downstream re- 
gion corresponds to x > 0. Hereafter, we will label with the subscript 0, 1, 2 
quantities calculated at xq, x = 0~ and x — + respectively. 

An implicit solution of Eq. [T]in the upstr eam region may be found by gener- 
alizing the approach used by iMalkovl (j 19971 ) . Eq. [1] can be spatially integrated 
in the upstream region from x to xq leading to: 



D(x,p) 



df 
dx 



u{x)f{x,p) = 



D(x,p) 



df 

dx 



1 



(2) 

The solution of the homogneneous equation associated to Eq. [5] still reads 
fh(x,p) = exp [ip(x,p)} with 



ip(x,p) 



D(x',p) 



(3) 
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from which the general solution follows as: 
/(*,p) = /i(p)e*p[#r,p)](l+«o [ dx' en)hr{1 '- l>)] 



D(x',p) 



uohip) 



+ Z(x',p) 



where 



Z(x,p) 



uofi(p) 



dx 



r du(x') 
dx' 



fi{p) — /(0,p) is the particle spectrum at the shock location, and 

4>esc(p) = - 



D(x, P ) d / x 



(4) 
(5) 

(G) 



is the flux of particles escaping from the shock across the surface at x = x 
(escape flux). This solution does not explicitly show that f(xo,p) = 0, although 
this condition has been clearly used in passing from Eq.[T]to Eq.[2] On the other 
hand, the condition /(xo) = directly leads to an interesting expression for the 
escape flux (f> esc (p), as soon as the transport equation is integrated between xq 
and the shock location: 



".,Ap) : { 1 - "n / dx ^L^^ P ^ z i x > P) 



exp[-ip(x,p)] 



(7) 

It is finally convenient to introduce the dimensionless functions K(x,p) and 
W{x,p) defined respectively as 



K{x,p) — uq dx 



D(x',p) 



-Z{x',p) 



W(x,p)=u ,V eXpH/;(a;>)] 



and 

D(x',p) 

so that the solution of the transport equation becomes: 

W{x,p) 



f{x,p) =/i(»exp [il>{x,p)\ ^l + K(x,p) 
and the escape flux: 

</w(p) = -wo AO) 



W (p) 

l+K (p) 
W (p) ' 



[1 + K (p)] 



(8) 



(9) 



(10) 



(11) 
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At this point we follow the procedure described e.g. in iBlasil (|2002l I2004T ) of 
integrating Eq. [T] across the shock and between xo and 0~ in order to derive an 
equation for flip), i.e. : 



Pdfi, 



hip) 



p du p 
3~dp 



+ 4>esc(p) ~ Ql(j>) 



(12) 



where we introduced the mean velocity effectively felt by a particle with mo 
mentum p in the upstream region: 

1 f° . du(x) 



u p {p) = ui 



dx- 



dx 



-f(x,p)- 



(13) 



Following (see iBlasi. Gabici fc Vannonil . 120051 ) we write the injection term as 

(14) 



Q(x,p) = Qi(p)5{x) = 7 ]™°2° 8(P ~ Pinj)6(x) 



where r\ is the fraction of particles crossing the shock and injected in the 
acceleration process, and jPj ra7 - is the injection momentum. As discussed by 



Blasi. Gabici fc Vannonil (|2005| ). we write the injection momentum as a multi- 
ple of the thermal momentum of particles downstream, pi n j — ^i n jPth,2- In the 
assumption that the thermal particles downstream have a Maxwellian spectrum, 
the fraction r\ is uniquely determined by the choice of £mj- 

We also introduce the normalized fluid velocity U — u/uq and the normalized 
escape flux $ esc — 4> eS c/ (uofi), as well as the compression ratios at the subshock 
Rsub — U1/U2 and between Xq and downstream R tot = u$/u2- The solution of 
Eq. [12] then reads: 

V n 3R tot I f p dp'3R tot [U p (p')-<S> esc {p')}\ 



A(p) = 



47rp?„j RtotU p (p) - 1 



exp ■ 



Pinj 



P > 



RtotUpip') - 1 



(15) 

It is easy to check that in the t est-particle limit K(x,p) = , U p {p) = 1 



and the standard solution (see e.g. 
in eqs. [IB [n]and [15] 



Caprioli. Blasi. Amato 



2009) is recovered 



3. An approximate solution 

In this section we use a heuristic argument to derive an approximate solution 
of the problem. An a posteriori comparison with the exact solution derived 



G 



above shows an excellent agreement in all cases considered. 

Let us consider the function K(x,p) defined in Eq. [5] At any given mo- 
mentum p, the distribution function can be regarded as approximately con- 
stant (f(x,p) ~ fi{p)) for x < x p ~ D(p)/uix p ) and exponentially suppressed 
for x > x p . Hence, for x <C x v we have Z(x,p) ~ U(x p ) — U(x) < 1 but 
also i/j(x,p) <C 1. This leads to K(x,p) ~ (uox/D(p))Z(x,p) <C 1. On the 
other hand, for x 3> x p , Z(x,p) — > so that both exp [ip(x,p)j K(x,p) and 
Kq(p) /Wo(p) tend to 0. For x ~ x p , the situation is less clear, but one can ex- 
pect that Z(x,p) <C 1 because a; — > x p and /i(p) starts feeling the exponential 
suppression. This suggests that we can neglect K(x,p) with respect to unity in 
Eq. I10[ although clearly this conclusion needs to be checked a posteriori against 
the exact solution. 

The following recipe is thus proposed as an approximation of the exact so- 
lution of the transport equation: 

w(x, P y 



f(x,p) = /i(p)exp 



4>esc(p) = 



D{x',p) 
uofiip) 



1 



W (p) 



(16) 



(17) 



W (p) 

This expression tends to the correct test-particle limit, as one can easily verify. 
A point that is worth highlighting is that the distribution function at the shock, 
fi(p), is sensible to the assumed spatial dependence only through the function 
Up, and is therefore weakly affected by whether the approximate or the exact 
solution is adopted. 

We want to stress an important point: in the case of boundary condition 
in momentum, namely when the maximum momentum is fixed, the procedure 
above would lead us to the functional form f(x, p) = fi (p) exp — dx' K , 
which is slightly different (and simpler) than any ansatze previously propose d 
in the literature (see for instance 



Malkov 



1997 



Blasi. Amato fc Caprioli . 



2007J). 



It turns out that the approximation found here, which is the simplest possible 
extrapolation of test-particle theory, gives a solution that is basically undistin- 
guishable from the exact one. In other words, both in the case of boundary 
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condition in momentum (fixed p ma x) or free escape boundary (fixed Xq), the 
best description of the spatial distribution of accelerated particles is provided 
by the simplest possibility which automatically satisfies all the relevant limits 
and the boundary condition at the shock. 

A full solution of the system of conservation equations for mass, momentum 
and energy, coupled with the diff usion-convection one, is obtaine d following the 



iterati ve procedure described by lAmato fc Blasil ([20051 ) and by ICaprioli et al 



(|2009t ) . when the gene ration of magnetic turbulence via streaming instab ility 
(see e.g. 



Skillina 



19751 ) and its dynamical feedback ([Caprioli et al 



2008) are 



taken into account. The only difference here is that there is no need to fix a 
maximum momentum by hand, since the distribution function gets intrinsically 
suppressed above a certain p max as a consequence of the escape at x = xq. 

The iterative method can be summarized as follows. Let us consider the 
momentum conservation equation, normalized to pou 2 ,: 

1 



U{x) + P c (x) + P w (x) +P g (x) = l + 



7 M 2 ' 

where we introduced the normalized pressure in cosmic rays: 

47T 



Pc{x) 



dpp v(p) f(x,p) 



(18) 



(19) 



the normalized pre ssure in magnetic tur bulence generated via resonant stream- 
ing instability (see 



Caprioli et al 



Pw(x) 



20091 Eq. 42): 
v A l- U(x) 2 



iu U{xf/ 2 ' 

and the normalized pressure of the background gas with adiabatic index 7: 

U(x)-~> 



(20) 



7 M 2 



(21) 



The last expression holds provided the heating in the precursor is purely adia- 
batic: its general ization to cases with s ome turbulent heating is however straight- 
forward (see e.g. ICaprioli et al.l . l2009l sec. 6). 

We start from a guess value for U% = R S ub/Rtot, which uniquely determines 
Pwi, Pgi and Pd via Eqs . l20l [2T1 and [TBI We notice that, once U\ is fixed, R su b 
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and Rtot can be worked out separately from the conservation equations in the 
precursor obtaining ( see ICaprioli et al.1 . 12003 . eq. 16) 

+1 = M$Rl ub 7 + l-^ sttb ( 7 -l) 

tot 2 l + P„i/^[l + ^U(2/7-l)] 

and in turn 



7? _ ^-V^-8^i(7 + l)(7-2) . 

n tot — ATT „ / , -K su b — U\ti to t \to) 

4(71-^,1(7 ~ 

where A = 2 7 (P gl + P wl ) + U l { 1 - 1). 

At this point, we start with a test-particle guess for f(x,p), normalized in 
order to account for the obtained P c i, and calculate P c (x) by using Eq. [T9l and 
then U (x) through Eq. [TU This updated velocity leads to a new P w (x) and 



hence to 5B{x) = y 8ttpqUqP w (x) which is used to update D(x,p). 

Now we can calculate a new f(x,p) as a function of the old distribution 
function and of the new U(x) and D(x,p), according to Eq. [Pol and Eq. [10] (or 
Eq. I16p. The new f(x,p) is again normalized to P c \ and the procedure above 
is iterated until convergence is reached, i.e. until f(x,p) and its normalization 
factor do not change between two successive steps. 

For an arbitrary value of Ui, however, the required normalization factor will 
be different from 1, thus the process is restarted with a different choice of U\ 
until no further normalization is needed. The distribution function calculated 
with the value of R su b /Rtot obtained in this way is by construction the solution 
of both transport and conservation equations. 

We consider here two cases: a test-particle-like one (inefficient acceleration) 
and a strongly modified one. In the first case, we choose £j„j = 4.3, corre- 
sponding to a fraction of injected particles r\ ~ 1.7 x 10~ 6 , while in the second 
case we choose ^i n j = 3.3, corresponding to rj ~ 1.2 x 10~ 3 . Moreover, in 
the inefficient case we assume Bohm diffusion in the background magnetic field 
Bq = 5/iG, while in the strongly modified case, having in mind the case of 
shocks in SNRs, we adopt a Bohm- like diffusion coefficient calculated in the 
magnetic field which i s self-generated t hroug h resonant streaming instability by 



accelerated particles (jAmato fc Blasi . 



2006J), namely D(x,p) = f^f^j. The 
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1Q -2 1Q 1Q 2 1Q 4 1Q 6 

P*= P/ m c 



Figure 1: Particle spectra at the shock and escape flux in the test- particle-like case multiplied 
by a factor 100 (§ = 4.3) and in a strongly modified case (£ = 3.3). Symbols correspond to the 
approximate solution given by Eq. 1161 while lines correspond to the exact solution (Eq. HOjl . 

other parameters are chosen as follows: the shock velocity is uq =5000 km/s, 
the free escape boundary is located at xo — 0.15 pc, and finally the background 
density (temperature) is po = 0.1m p cm -2 (To = 10 5 K), corresponding to a 
sonic Mach number Mq ~ 135. Again, these choices are inspired by the values 
expected in SNRs. 

In Fig. [T] we plot the particle spectrum at the shock and the escape flux in a 
test-particle-like case (multiplied by a factor 100, lower curves) and in a strongly 
modified one (upper curves). In Fig. [2j instead, we show the hydrodynamical 
quantities in the upstream region (top panel), and the distribution function at 
some given upstream positions, namely at x/xq = 10~ 7 , 10~ 4 , 10~ 2 , 0.5 (bottom 
panel), all referred to the modified case. All the curves in the Figs. [1] and [2] 
refer to the exact solution. The approximate solution given by Eq. [T5] leads to 
the results shown with symbols in both figures. One can easily realize that the 
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Figure 2: Upstream hydrodynamical quantities (top) and cosmic ray distribution function 
at x* = x/x = 1CT 7 , 1(T 4 , l(r 2 ,0.5 (bottom), in the £ = 3.3 case. In any panel symbols 
correspond to the approximate solution and lines correspond to the exact solution. 
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agreement between the results obtained with the exact and the approximate 
solution is excellent, beyond any expectation. Moreover, the case of inefficient 
acceleration reduces exactly to the test-particle case, as shown by the lower 
curves in Fig. [T] 

The efficient case shows the typical features of cosmic ray modified shocks, 
with a c oncavity in the spectrum induced by the precursor in the upstream fluid 



(see e.g. iMalkov fc Drurvl . I2001L for a comprehensive review). In the present 
case the total compression coefficient is Rtot — 10.6 as a result of the pressure 
in cosmic rays (about 66 per cent of the bulk pressure at the shock) and of the 
dynamical backreaction of the amplified magnetic field, since upstream the mag- 
netic pressure domin ates over the gas pressure {P w i/Pgi — 45), as described in 



Caprioli et al.l (|2008f) . In this case the energy carried away by escaping particles 
represents about 37 per cent of the bulk energy flux. 

It is interesting to notice that in the case of efficient acceleration, despite the 
fact that the magnetic field amplification induced by accelerated particles at the 
shock is of order SB/Bq ~ 20, the maximum momentum which is implied by 
the free escape boundary condition at x — xq is only a factor ~ 2 higher than 
in the inefficient case (we recall that in this latter case the diffusion coefficient 
is Bohm-like in the background magnetic field -Bo)- This apparently counter 
intuitive result is in fact simple to understand: due to the dynamical reaction 
of the accelerated particles, the effective fluid velocity felt by particles in the 
precursor is U\ ~ 0.3, which implies a slower acceleration rate and lower maxi- 
mum momentum; moreover the fact that in the efficient scenario the magnetic 
field is self-generated implies that most of it is concentrated around the shock, 
while the (turbulent) magnetic field responsible for particle diffusion close to 
xq is in fact much smaller than Bq (SB turns out to be smaller than Bq for 
x < 0.5iEo). As a consequence of these two facts, the maximum momentum does 
increase in the modified case with respect to the inefficient one, but less than 
the naive expectation would suggest: in fact, as far as our investigation of the 
parameter space has gone so far, the maximum momentum does scale linearly 
with the position of the free escape boundary xq but not with the strength of 
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the amplified magnetic field. 



4. Discussion and Conclusions 

Here we discussed the first semi-analytical exact solution of the problem 
of particle acceleration in non-relativistic shocks with a free escape boundary, 
when non-linear effects induced by the dynamical reaction of accelerated par- 
ticles and by the amplification and dynamical feedback of the magnetic field 
are taken into account. In addition to the exact solution, which is rather cum- 
bersome to implement in numerical calculations, we also proposed a simple 
but excellent approximation to the exact solution. This approximate solution 
catches all the main Physics ingredients of the problem and is computationally 
very convenient. We checked this approximate solution versus the exact solution 
and the agreement, both in terms of the spectrum of accelerated particles at 
the shock and in terms of the spatial distribution of accelerated particles in the 
precursor, is excellent. As a consequence, also the shock structure in terms of 
spatial dependence of the hydrodynamical quantities and of the self-generated 
magnetic field is perfectly reproduced. The escape fluxes and spectra are also 
in stunning agreement. 

The ability at providing not only the spectrum of accelerated particles, but 
also the spectrum of particles escaping through the free escape boundary lo- 
cated at a position xq upstream, is exactly what makes the solutions presented 
here (both the exact one and the approximate one) especially valuable. In a re- 
alistic situation, such as the expanding shock front associated with a supernova 
remnant, the existence of a free escape boundary leads to a maximum momen- 
tum of the accelerated particles which depends on time and in general decreases 
with time during the Sedov- Taylor phase, if the magne tic field is generated by 



the a ccelerated particles through streaming instability (jCaprioli. Blasi. Amato 



20091 ). It follows that the convolution in time of the instantaneous escape flux 
leads to the formation of a complex spectrum which is no longer peaked around 
a specific momentum. 
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From the physical point of view the main uncertainty related to this type of 
calculation is 1) in the nature of the self-generated waves and their interaction 
with accelerated particles and 2) in the determination of the location of the free 
escape boundary based on first principles. These two issues, clearly related to 
each other, are not easily solvable at the present time and a phenomenological 
approach is the only one we can afford to adopt. 

From the mathematical point of view, the solution presented here is an 
important step forward in the description of the process of particle acceleration 
in astrophysical shocks, especially in SNRs. The limitation that remains is 
that the solution assumes that the system is able to reach a quasi-stationary 
configuration at any given time. Despite this limitation these methods are 
of the greatest importance in order to have an appropriate description of the 
acceleration process in complex astrophysical objects such as SNRs. Other 
methods, all numerical in nature, have in fact typical running times that range 
between several hours and several days for a given set of parameters, compared 
with O(minutes) required by the semi-analytical approach presented here. 
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